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INTRODUCTION TO THE APPLICATION OF 
YON ZEIPEL'S METHOD 

by 

William J. Wickless, Jr. 

INTRODUCTION 

To understand the application of Von Zeipel's method in the case of 
a satellite orbiting under the influence of a gravitational field plus small 
perturbating forces, it is first necessary to summarize the development 
of several concepts of classical mechanics. Therefore, the first portion 
of this paper will be a brief resume of material which may be found in 
greater detail and slightly modified form in Chapters 7-8 of Classical 
Dynamics of Particles and Systems by Jerry B. Marion, Dept, of Physics, 
University of Maryland. A good knowledge of ordinary differential and 
integral calculus is assumed. 

I. AN INTRODUCTION TO THE CALCULUS OF VARIATIONS 

The basic problem of the calculus of variations is determining the 
unknown function y(x) describing a path between two fixed end points 
we call y(Xj ) and y( x 2 ^ such that 

P 

I = f [y(x), y'OO, X J dx 


1 


takes on a maximum or minimum, f is a given function of the functions 
y(x), y'(x) and the independent variable x and the limits of integration 

are fixed. That is, if we wish to minimize I (say), we wish to find a 
function y(x) such that if y x (x) is any other continuous function such 
that YjfXj) = y(Xj), y x (x 2 ) = y(x 2 ) - any other path between y(Xj) and 
y(x 2 ) - then: 

r 2 r 2 

f [y(x), y'(x), x] dx < f [y^x), y'(x), x] dx 



Xj x 2 assumed to be differen- 

tiable to any needed 
order) 

Two possible paths between y(x x ) and y(x 2 ) are sketched above. 

We will consider families of functions, giving possible paths between 
two fixed points y(x^ and y(x 2 ), indexed by a parameter a, a running over 
a suitable segment of the real line. That is, we will consider families 
of the form 


xe [x Jf x 2 ] 

{y(a. x)} 

ae [a, b] with 


2 



"y (a, x x ) = y(Xj) = constant' 


y(a, x ) = y(x ) - cons t an t 


for all 
ae [a, b] , 


Then the integral I becomes a function of the parameter a 


For example 


where 


i 


1(a) = f [y(a, x), y' (a, x), x] dx . 


2tt 

J * 


1(a) = | x (a sin x) dx 

■'0 


ae [1, 2] 


y( a - x) i a sin x 

f [x, y (a, x), y' (a, x)] = xy (a, x) 

Xj =0, X 2 = 277 
y ( x j ) = o y(x 2 ) = 0 

NOTE: a sin (0) = 0 ; a sin 2 tt = 0 for all ae [1, 2] ; {a sin x} 

sents the family of paths between the pts (0, 0) = y(Xj), (277, 
shown at top of following page. 


repre- 

) = y( x 2 ) 


3 


■»• y(3/2, x) = 3/2 sin x 


y(XjV 





1 





y(x 9 ) 


Now if the integral so written is to have an extremum along some 
path y(a Q , x) 

a. Q e [a, b] 

then 


31(a) 

3a 


a=CL o 


= 0. 


This means 


3a 


x 2 


f [y(a, x), y' (a, y), x] dx 


= 0. 


x i 


a=a 


0 


since x is independent of a, the partial may be taken under the integral 
sign: 
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= 0 . 


x 2 



f 

"Bf By Bf By* 

dx 

J 

By Ba By' Ba 


X 1 




Now the above equation may be integrated by parts (see (1), Sec. 7-3) to 
obtain 


x 2 



f 

Bf d _B_f 

dx 

J 

By dx By' 

Ba 

X 1 




As f is assumed to be an arbitrary function of y(a, x) y' (a, x) x , 
the integrand itself must vanish: 

Bf d Bf By _ q 

_3y dx By'J Ba 


for all xe (Xj, x 2 J. 

Also since {y(a, x)} is an arbitrary family of paths, in general 


By (a, x) 
Ba 


* 0 

a " a 0 


along the curve 


y = y ( a 0 . x ) 


X€ [Xj, x 2 ] 
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Now if x(t) is the path followed by the particle, the Euler - Lagrange 
equation must be satisfied, i.e. 


BL d_ BL 

Bx(t) dt 3x(t) 

Now 


BL 

Bx(t) 


Kx(t); 


BL 

Bx(t) 


= m x(t) 


- Kx(t) = ~ [m x(t)] 
d t 

or 

x(t) + w2 x(t) = 0 with w Q = /K/m . 

This differential equation may be solved (see any elementary text) 
to obtain 


x(t) = A sin (w Q t + S) 

where A + S are constants of integration to be fitted to the initial data - 
in this case A = h , 8 = n/2 . 

Now an important thing to notice is that the Lagrangian is a scalar 
function (being the difference of kinetic and potential energies - both 
scalar functions). Thus we can, if desired, transform one co-ordinate 
variable x(t) to a new variable q(t) under some co-ordinate transfor- 
mation and we have 


8 



L(x(t), x(t), t) =L(q(f), q(t), t) . 

(This is actually the defining property of a scalar function). 
The Euler -Lagrange equation is still valid: 

3L _ _d_ c»L _ Q 
3q(t) " dt Bq(t) 


Now as one final generalisation, it is clear the original analysis 
could have been carried out to treat the case of a particle with more 
than one degree of freedom. The Lagrangian then would have been a 
function of several co-ordinate variables, their derivatives, and pos 

sibly the time, i-e. 


L = L(Xj(t) 


■x n (0. X x (t-) x n (t), t) 


We then could proceed in a similar fashion 
instead of the single equation 


[see (1), Sec. 7-7] to obtain, 


BL d BL _ Q 
3x dt 3x 


the system of equations: 


BL_ _ _d_ BL. _ o 
Bx j ” dt Bxj 


i = L 


2 


N. 
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Now as this Lagrangian is still a scalar function, we may make if 

desired, a co-ordinate transformation 0 from the position variables 

^ x r/ t )} to some new set of variables (q^t) . . . . q 
and we still have: 


Lfaj 

And also 



O =L(x, 


t). 


3L d 3L 

3qj dtBqj - ^' j - 1 . . . n (a) 


as a 


It is important to note that the new variables may be chosen 

— y_ N functions q, C t> - q n (t) at all as long as the Lagrangian 
may be expressed as L = L(q,(t). ■ ■ qJt ), q l( t) . . t) . In 

a real problem, a set of variables • -q^t) are usually chosen 

so that L may be easily expressed as Uq, • • • q„ q, . . . ^ t) . Then 

the system (a) is integrated, giving qj vqj (t , constants of integration) 
j = 1 ’ • ’ n (b) . 


Finally the transformation 0 is inverted and 0-i is applied to (b) 
yielding 


xj = xj(t) j = 1 N. 


The vector 
the particle unde 


- (Xjft) • • • • x n ( t)) is then the 
r consideration. 


position vector 


of 
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= 0. 


x 2 



f 

Bf By Bf By’ 

dx 

J 

By Ba + By' Ba 


X 1 




Now the above equation may be integrated by parts (see (1), Sec. 7-3) to 
obtain 


Bf 

By 


_d_ _3f 
dx By' 


^ dx 

da 


As f is assumed to be an arbitrary function of y(a, x) y' (a, x) x , 
the integrand itself must vanish: 


Bf d Bf 
By dx By' 


By 

da 


= 0 


for all xe [Xj, x 2 ] . 

Also since (y(a, x)} is an arbitrary family of paths, in general 

£ 0 along the curve 

a ~ a o 


By(a, x) 
Ba 


y = y (%’ x ) 

xe [x r x 2 ] 


therefore 


Bf _d_ Bf 

^y( a 0 - x) dx 3y'(a 0 , x) 


for all xe [x x , x ] . 

This is known as the Euler - Lagrange equation. 

II. APPLICATION OF THE EU LER -LAGRANGE EQUATION 
TO PHYSICAL SYSTEMS 

It is known that if a particle is allowed to move between two fixed 
pts in the plane x(tj') and x(t 2 ) in the time interval [t x , t ] in any of 
the possible paths given by the family of curves {x(a, t)^ te [t Jt t ] 

ae [a, b] 

it will actually move in the path x = x(a Q , t) for which 


(T - U) dt is a relative extreme. 

(Note that this is a local 
condition. ) 


(see ( 1) , Sec. 8.3) 


Here T is the kinetic energy of the particle at the point x(a, t) - 
T - 1/2 m x 2 (a, t) — and U is the potential energy at the same pt 
U = U(x(a, t)). 


6 



along the path of the particle's motion x(t) = x(a Ql t) , the Euler - 
Lagrange equation must be satisfied by the quantity L = (T - U). L is 
called the Lagrangian of the system. 

Thus we have 

BL _d BL _ 

"dt 2x'(t) ~ 

If the Lagrangian of a system is known, the above equation may be 
integrated giving the path of the particles motion x = x(t) . 

EXAMPLE: Determination of the Equation of Motion for a 

Harmonic Oscillator 

Let a particle of mass m be attached to a spring with spring constant 
K and relaxed length 0. Choose a co-ordinate system such that the point 
of equilibrium of the spring is at x = 0, Let the mass be pulled to a 
distance h above the x = 0 level and released at a time t = 0, Find x(t) 
(the bar is hereafter dropped for convenience) 


x = 0 

Now for any time t, we have T = 1/2 m x 2 (t) and U = 1/2 K x 2 (t) . 
so L = T- U=l/2mx 2 (t) - 1/2 K x 2 (t) . 
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Now if x(t) is the path followed by the particle, the Euler -Lagrange 


equation must be satisfied, i.e. 


BL _d_ 3L 

3x(t) dt 3x(t) 


Now 


3L 

3x(t) 


Kx(t); = in x(t) 

3x(t) 


" " Kx(t) = [m x(t)] 
d t 

or 

x(t) + w2 X(t) n 0 with w Q = /K/m . 

This differential equation may be solved (see any elementary text) 
to obtain 


x(t) = A sin (w Q t + 5) 

where A + § are constants of integration to be fitted to the initial data — 
in this case A = h , S = 77 / 2 . 

Now an important thing to notice is that the Lagrangian is a scalar 
function (being the difference of kinetic and potential energies - both 
scalar functions). Thus we can, if desired, transform one co-ordinate 
variable x(t) to a new variable q(t) under some co-ordinate transfor- 
mation <$> and we have 


8 



L(x(t), x(t), t) = L (q( t), q(t), t) . 


(This is actually the defining property of a scalar function). 

The Euler - Lagrange equation is still valid: 

BL _ _d_ BL = 

3q(t) dt Bq(t) 

Now as one final generalization, it is clear the original analysis 
could have been carried out to treat the case of a particle with more 
than one degree of freedom. The Lagrangian then would have been a 
function of several co-ordinate variables, their derivatives, and pos- 
sibly the time, i.e. 


L = L(Xj(t) * • ■ -x n (t), x x (t) • • • • x n (t), t) . 

We then could proceed in a similar fashion [see (1), Sec. 7-7] to obtain, 
instead of the single equation 


BL BL 
Bx d t Bx 


the system of equations: 


BL _d_ BL 

Bxj dt Bxj 


j = 1, 2 N. 
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Now as this Lagrangian is still a scalar function, we may make if 


desired, a co-ordinate transformation^ from the position variables 
|xj(t)*** x n (t)j to some new set of variables {qj(t) .... q n (t)j- 
and we still have: 

L(qj • ’ ’ ’ • q n . q, q n - t) = L(x x • • • x n , ^ • x n , t). 


And also 


BL _d_ oh 

Bqj dt Bqj 


j = 1- • • n (a) 


It is important to note that the new variables may be chosen as a 
set of any N functions q x (t) • • • • q n (t) at all as long as the Lagrangian 
may be expressed as L = Lfqj ft) • • • q n (t), q t (t) • • • q n (t) t) . In 
a real problem, a set of variables qjft) 1 ■ ’q n (t) are usually chosen 
so that L may be easily expressed as Lfqj • • • q n cjj • ■ • q n t). Then 
the system (a) is integrated, giving qj =qj(t, constants of integration) 
j = 1 • ■ -n(b). 

Finally the transformation is inverted and c£ -1 is applied to (b), 
yielding 

xj = xj (t) j = 1 • ■ ' * N. 

The vector rft) = (x ft) ■ • • ■ x n ( t) ) is then the position vector of 
the particle under consideration. 


10 



III. CONSTRUCTION OF THE HAMILTONIAN; THE 


CANONICAL EQUATIONS OF MOTION 


Now let us consider a system in which (1) the Lagrangian is not an 
explicit function of the time , i.e. L = L(q (t) • * q n ( t), q t (t) * • * q n (t)) 
and (2) the q. ! s have been obtained from the x.'s by some transforma- 
tion cjj that does not explicitly involve the time. 

Let us define N new quantities: 


P 

) 


3L 

3<i i 


j =1,2 N. 


(I) 


P is called the conjugate momentum associated with q . 

Using this definition, the Euler - Lagrange equations for our system 
may be rewritten in the form 


FV = ^ j = 1, 2 • • • - N. (II) 

Now let us define a new quantity 

h 'Z] p ^, - l - 

j 

H is called the Hamiltonian of the system, and under assumptions 
(1) and (2) above, it can be shown that H is equal to the total energy of 
the particle under consideration. H = T + U* (See (1), Sec. 8-8). 
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Now if H is considered to be a function of q^.Pj j - 1 • * 'N, then 


dH 


= T is.dc. t |S-dP. 

L—i Bq . BP. J 


also since 


H 


■ L p a - 


(A) 


dH 


T P. dq. +q dP. dq dq 

/ i J J J 1 dqj 1 
j 

using (I) and (II), the above becomes 


dH 


y [P. dq. - P j dq.] + q j dP j - P j 


. dq. 
j J 



(B) 


Now identifying the coefficients of dP. & dq. in (A) and (B) we have 


P _j^L - BH 

i Bq . ^ BP. 

J 


j = 1, 2 • 


N. 


These are called the cyclic equations of motion. As in the case of 
the Euler -Lagrange equations, they may be integrated to obtain the 
equations of motion of the particle under consideration. 


12 



We will now consider our problem — determination of the equations of 
motion for a satellite orbiting under the influence of a gravitational field 
plus small perturbations (e.g. the effects of the earth’s oblateness.) The 
following (sec. IV) is a summary of a portion of the material found in 
the paper "Notes on Von Zeipel’s Method" by Giorgio E. O. Giacaglia, 
published at the Goddard Space Flight Center. 

IV. VON ZEIPEL’S METHOD 

In the case we are considering, the Hamiltonian may be expressed 
as a function of six variables (L, G, H, Z, g, h) where L, G, H are the 

conjugate momenta as sociated with £ , g, h respectively. These six 

variables are called the Delaunay variables. 

The equations to be integrated for our problem then are: 

} BH . BH / BH 
'i = — g = — h = — 

BL 3G BH 


_^ BH £ = BH £ = BH 

B'fc 3g Bh 

here we have just taken the cyclic equations of motion and set L = P 1 , 

G = P 2 , H = P 3 , l = q t , g = q 2 , h = q 3 . 

In many cases the negative of the Hamiltonian F = - H is introduced; 
the cyclic equations of motion for our system then become: 


13 


BF . 3F ^ _ BF 
BL g - " BG BH 


L 


BF 

B^ 


G = 


BF 

Bg 


H = 


BF 

Bh 


Now in solving the above system, it would greatly simplify our 
problem if we could transform our variables (L, G, H, I, g, h) to a new 
set of variables (L’ , G' , H ' , V , g' , h' ) in which the new function F could 
be expressed in a form not explicitly containing one of the new variables 
(-£' say), i.e. 


F' = F'(L', G', H' _ g\ h') 

then 



V 


cons tant . 


Here we must remember that, for the cyclic equations of motion to 
be valid, the new variables must satisfy 


L' 


BL 

W 


G' 


BL 

Bg' 


H' 


BL 

Bh' 


i.e. L', G' , H' must be the conjugate momenta with respect to <t , g' , 
h' . (This assumption was used to derive the cyclic equations of motion). 

If the above equations are satisfied, then ( L' , G' , H ', ■i , g' , h' ) are 
said to be a canonical set of variables, and the transformation connecting 
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the primed and unprimed variables is called a canonical transformation 


[assuming F(L, G, H, l, g, h) = F'(L' } G', H\ - g', h' ) for all times t]. 

Now finding the desired canonical transformation connecting ( L, 

G, H.I.g, h) to a new set of variables (L', G', H', -V , g', h' ) for which 
BF'/B-t' = 0 is equivalent to finding a function S = S(^,g, h, L', G H'), 
called a generating function, satisfying: 

L = — G = — H = — 

U Bg Bh 

r _ ds_ , _ BS h , _ BS 
BL' g " BG' " BH' ’ 


This function S will determine an implicit transformation connecting 
(L, G, H, -t , g, h) and (L‘, G\ H', -t', g', h' ) and it is known that the new 

set of variables will be canonical. (See Planetary Theory by Brown and 
Shook, sec. 5-3.) 

We now shall give explicit formulas for determining S and F' . Here 
we first employ the assumption that our satellite is orbiting in a gravi- 
tational field, subject to small perturbing forces. We assume the 
function F (L, G , H , ^ , g , h ) can be expanded in a series 


F 


= F o + 
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where F Q is the negative of the Hamiltonian calculated in the case of a 
body orbiting under the influence of gravity alone and X is a small 
parameter . (For the calculation of F Q = fj? / 2L 2 where /x — K M , with K 
the Gaussian constant and M the mass of the orbited body; see any 
intermediate textbook in celestial mechanics.) 

We assume similarly: 

00 

= f ; ♦ 2Z F i x ‘ 

J = 1 

CD 

s = s » + 22 s - xi 
j = l 

where S 0 is the generating function determining the identity transforma- 
tion, S 0 = VI + G'g + H'h. (In the Keplerian case the original set of 
variables are the ones desired; no transformation is necessary.) 

Now to any order in k , F (L, G, H, 't, g , h ) = F' (L', G r , H' - g , h ) . 
(Remember F' is not to be an explicit function of ^ .) 

Look at this equation up to second order in i.e. 

F 0 + xf 1 + x 2 f 2 =f; + xf; + \ 2 f; (A) 


now 
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bs y> j 

L = B^ = U + L> 

j=l 



j=l 


with similar expressions for G, H, h • 

Substituting the above into (A), we have 



| 

1 

1 

■ 
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G' + 





G', H ' — , 



BS, 

BG 1 



Now expanding the left side of (B) in a Taylor series around L, G, 
H and Fj and F' on the right hand side in a Taylor series around g' , h' 
and keeping only terms up to and including the second order in X, we 
have 


F 0 ( L ') + 



BS, 

BS, 

1 

X + — X 2 

_ B-C 

d'i 


+ \ 2 F x (L', G', H', l, g, h) +\ 2 


pFj 

B Sl 

B Fl 

BS] 

BF X 

BSl 

_BL 


BG 


+ BH 

Bh_ 


+ ^ 2 F, [L', G', H', l, g, 


h] = f;(l', g', h' - g \ h'> 
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+ \F| (L', G', H' - g, h) + X 2 


'BF/ 

BS X 

BF, 

BS, 

L*v 

BG' 

+ Bh' 

BH'. 



[L', G', H' - g, h] 


(C) 


Now equating the coefficients of like powers of \ in (C) we obtain 


F!(L’,G',H'-g',h')=F 0 (L") = i r3 

2L 


(D) 


Fj (L' , G', g, h) + 


BF Q BS 

1l 


Hi = F; (L' , G' , H' - g, h) 

u 1 


(E) 


BF 0 BS 2 BFj BS, BF, BS, BF, BS, 

BL p^. BL BG BG' BH Bh 


bf; bs, bf; bs, 

+ F 2 (L\ G'.H'H, g, h) = — — , + 


+ F' (L'.G', H' - g, h). 


(F) 


Now F; is determined by (D) so (E) may be solved by putting F| - 
part of F, independent of l. This means (BF/BL^BS/B^ = - (? t - F j ) 


.19 


(part of Fj dependent on Z). Now F' and Sj are determined with the 
desired properties, and therefore (F) may be solved for S 2 by putting 
F 2 = part of F 2 independent of 


This determines F 2 and S 2 with the required properties. This 
process may be continued to determine F n ' and S n to any desired order. 


' j 
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